PLN et EMtree pour REMMOA, sans spatial

library(EMtree)

## Warning: replacing previous import 'dplyr::combine' by 'gridExtra::combine'
## when loading 'EMtree'

library(tidyverse)
library(grid)
library(gridExtra)
library(tidygraph)
library(ggraph)
library(ggbeeswarm)

Chargement et mise en forme des données

• Lignes nulles
• Format des covariables
• Y_corr = comptages moyens par détection

setwd(WorkDir <- "/Users/raphaellemomal/these/Data_Oak_remmoa/REMMOA")
ShapeDir <- paste(WorkDir, "shape", sep = "/")
OutDir <- paste(WorkDir, "output", sep = "/")
theme_set(theme_light())
load(paste(OutDir, "20180620_PLN_REMMOANC.RData", sep = "/"))
colnames(X)[7:8]<-c("Xcart","Ycart")
null_index=which(rowSums(Y)==0)
Y=as.matrix(Y[-null_index,])
X=as.matrix(X[-null_index,])
N=as.matrix(N[-null_index,])

# correction des comptages par le nombre de détection
Y_corr=as.matrix(Y/N)
Y_corr[!is.finite(Y_corr)]<-0

Y %>%as_tibble() %>%  gather(species, count) %>%filter(count!=0) %>%  ggplot(aes(reorder(species,count), count, color=species))+geom_quasirandom()+labs(x="", title="Non nul counts")+ coord_flip()+
  theme(axis.text.x = element_text( angle = 45,hjust=1))+guides(color=FALSE)

# covariables numériques
X=data.frame(X)
X[,c(6:10,14:20)]<-apply(X[,c(6:10,14:20)], 2, function(x) as.numeric(x))

Xquanti=X[,c(6,9:10,14:20)] # sans Xcart et Ycart qui sont des transformations de longi et lati
plot(data.frame(Xquanti), pch=20, cex=0.5)

Fonctions utiles :

• Ajuster PLN avec des covariables de vec
• afficher le réseau après rééchantillonnage

get_model<-function(data, vec){
  t1<-Sys.time()
  string<-paste(deparse(substitute(data)), paste(vec, collapse=" + "), sep=" ~ ")
  formula<-as.formula(string)
  mat = as.matrix(lm(formula, x=T)$x)
  model<-PLN(data ~ -1+mat)
  t2<-Sys.time()
  print(difftime(t2,t1))
  return(model)
}
get_network<-function(resample_output,Pt,title, data=Y, f=0.8){
  df<-freq_selec(resample_output$Pmat,Pt)
  graph<-draw_network(df,title, pal="dodgerblue3", 
                              nodes_label=colnames(data), layout="nicely", curv=0.05,
                      filter_deg=TRUE)
  return(graph)
}

Modèles PLN

Plusieurs modèles PLN ont été ajustés et stockés dans models.RData.

• m0 : ~ 1
• m1 : ~ 1 + SEA_STATE
• m2 : ~ 1 + Depth
• m3 : ~ 1 + SUBJECTIVE
• m4 : ~ 1 + Depth + SUBJECTIVE
• m5 : ~ 1 + SEA_STATE + SUBJECTIVE
• m6 : ~ 1 + SEA_STATE + Depth
• m7 : ~ 1 + SEA_STATE + Depth + SUBJECTIVE
• m8 : ~ 1 + longitude + latitude
• m9 : ~ 1 + longitude + latitude + Depth
• m10 : ~ 1 + longitude + latitude + SEA_STATE + Depth + SUBJECTIVE

Les modèles sont ajustés pour les données brutes Y (pas Y_corr). Le temps d’ajustement varie entre environ 30s et plus de deux minutes pour m8 et m9.

• Dans PLN on veut des grandes valeurs pour BIC, ICL et loglik
• les modèles 2, 4, et 6 comprennent “depth”, qui semble essentielle
• les coordonnées spatiales seules (modèle 8) montrent de belles performances
• cependant le temps d’ajustement est très conséquent (2.6min contre 28 à 56s pour les autres)
• le gain en ICL de 10 par rapport à 7 est négligeable
• depth+sea_state est meilleur que depth+spatial

Comparons les réseaux obtenus avec m0, m2 et m7 (pas d’ajustement spatial donc).

Les réseaux

Totaux

Les espèces surlignées en jaune dans les réseaux sont celles avec un score de betweenness important. La matrice d’offsets est la multiplication de la longueur du segment d’un site et de la largeur effective de l’espece.

• R_brut : m0 et sans matrice d’effort. 5min, alpha mean = 0.6959584
• R_offset : m0 avec matrice d’efforts. 5.3min alpha mean = 0.710248
• R_depth : m2 et effort. 5.9min alpha mean = 0.7966224
• R_full : m7 et effort. 15.3min alpha mean = 0.9731524

load("results/ReseauxEntiers.RData")
naif=get_network(R_brut,"naif",Pt=NULL)
offset=get_network(R_offset,"Offset",Pt=NULL)
depth<-get_network(R_depth,"~ Depth",Pt=NULL)
depth_seastate<-get_network(R_depth_seastate,"~ Depth+sea_state",Pt=NULL)


grid.arrange(naif,offset,nrow=1, ncol=2)

grid.arrange(depth,depth_seastate,nrow=1,ncol=3)

Comparaison nord / sud

indicesNorth=which(X$Ycart > 7.76e6 & X$Xcart > 2e5 & X$Xcart<6e5)
indicesSouth=which(X$Ycart < 7.55e6 & X$Ycart > 7.4e6 & X$Xcart > 5.5e5 & X$Xcart<9.5e5)

YNorth=Y[indicesNorth,]
null_colnorth=which(colSums(YNorth)==0)
YNorth=Y[indicesNorth,-null_colnorth]
YSouth=Y[indicesSouth,]
null_colsouth=which(colSums(YSouth)==0)
YSouth=Y[indicesSouth,-null_colsouth]

load("results/ReseauxNorth.RData")
graph0<-get_network(R_North,"Offset", data=YNorth,Pt=NULL)
graph_effort<-get_network(R_north_depth,"~ Depth", data=YNorth,Pt=NULL)
graph_depth<-get_network(R_north_full,"~ Depth+sea_state", data=YNorth,Pt=NULL)

grid.arrange(graph0,graph_effort,graph_depth,nrow=1, ncol=3, top="Détections au Nord")

load("results/ReseauxSouth.RData")
graph0<-get_network(R_South,"Offset",data=YSouth,Pt=NULL)
graph_effort<-get_network(R_south_depth,"~ Depth",data=YSouth,Pt=NULL)
graph_depth<-get_network(R_south_full,"~ Depth+sea_state",data=YSouth,Pt=NULL)

grid.arrange(graph0,graph_effort,graph_depth,nrow=1, ncol=3, top="Détections au Sud")

Comparaison par région

Le plan d’échantillonnage en Nouvelle-Calédonie a été séparée en 2 zones : Océanique et talus continental, ces deux zones sont elles-mêmes séparées en trois zones de transects.

load("results/ReseauxRegion.RData")
load("results/DataRegion.RData")
Graphs=lapply(seq_along(ReseauxRegion), function(x){
  data=DataRegion[[x]]$Y
  get_network(ReseauxRegion[[x]], title=strsplit(names(ReseauxRegion)[x], split="_")[[1]][3], data=data,Pt=NULL)
})
grid.arrange(Graphs[[1]],Graphs[[2]],Graphs[[3]],ncol=3, nrow=1)

grid.arrange(Graphs[[4]],Graphs[[5]],Graphs[[6]],ncol=3, nrow=1)